#include <life/options.hpp>
#include <life/lifealg/backend.hpp>
#include <life/lifediscr/functionspace.hpp>
#include <life/lifefilters/gmsh.hpp>
#include <life/lifefilters/exporter.hpp>
#include <life/lifevf/vf.hpp>
#include <life/lifediscr/operatorinterpolation.hpp>
#include <life/lifefilters/geotool.hpp>

#include <iostream>

using namespace Life;
using namespace Life::vf;

typedef Application Application_type;
typedef boost::shared_ptr<Application_type> Application_ptrtype;

namespace opinterp
{
  
  template<uint Dim,uint OrderGeo>
  boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
  createMesh( std::string __name,std::string __str)
  {
    typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
    mesh_ptrtype mesh( new mesh_type );

    Gmsh gmsh;
    gmsh.setOrder(OrderGeo);

    std::string fname = gmsh.generate( __name, __str  );

    ImporterGmsh<mesh_type> import( fname );

    float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
    if (gmsh_version==2.0) import.setVersion( "2.0" );
    else if (gmsh_version==2.1) import.setVersion( "2.1" );

    mesh->accept( import );
    mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
    mesh->updateForUse();

    return mesh;
  }

  template<typename mesh_type,typename element_type>
  void
  exportResults( Application_ptrtype & test_app, element_type & U, std::string __str )
  {

    typedef Exporter<mesh_type> export_type;
    typedef boost::shared_ptr<export_type> export_ptrtype;
    export_ptrtype exporterLoc( Exporter<mesh_type>::New( test_app->vm(), test_app->about().appName()+__str ) );

    Log() << "exportResults starts\n";
    exporterLoc->step(0)->setMesh( U.functionSpace()->mesh() );
    exporterLoc->step(0)->add( __str, U );
    exporterLoc->save();
    Log() << "exportResults finish\n";

} // Interpolate::export


  template<uint Dim,uint OrderChamp,uint OrderGeo>
  void
  run(Application_ptrtype & test_app)
  {

    typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

    double meshSize1 = test_app->vm()["hsize1"].template as<double>();
    double meshSize2 = test_app->vm()["hsize2"].template as<double>();
    
    GeoTool::Node x21(-1,-1);
    GeoTool::Node x22(1,1);
    GeoTool::Rectangle R(meshSize1,"RIGHT",x21,x22);
    R.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
    R.setMarker(_type="surface",_name="Omega3",_marker1=true);

    std::string __name = "ola";
    mesh_ptrtype mesh = createMesh<Dim,OrderGeo>(__name,R.geoStr());

    GeoTool::Rectangle R2(meshSize2,"geo2",x21,x22);
    R2.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
    R2.setMarker(_type="surface",_name="Omega3",_marker1=true);

    std::string __name2 = "ola2";
    mesh_ptrtype mesh2 = createMesh<Dim,OrderGeo>(__name2,R2.geoStr());

    /*_________________________________________________*/

    typedef bases<Lagrange<OrderChamp,Scalar,PointSetFekete> > basis_type;
    typedef FunctionSpace<mesh_type, basis_type> space_type;
    typedef boost::shared_ptr<space_type> space_ptrtype;
    typedef typename space_type::element_type element_type;

    typedef bases<Lagrange<OrderChamp,Vectorial,PointSetFekete> > basis_v_type;
    typedef FunctionSpace<mesh_type, basis_v_type> space_v_type;
    typedef boost::shared_ptr<space_v_type> space_v_ptrtype;
    typedef typename space_v_type::element_type element_v_type;

    
    space_ptrtype Xh = space_type::New( mesh );
    element_type u( Xh, "u" );
    u = project( Xh, elements(mesh), Px()*exp(-2*Py()) );

    space_v_ptrtype Xhv = space_v_type::New( mesh );
    element_v_type uv( Xhv, "uv" );
    uv = project( Xhv, elements(mesh), vec(Px()*exp(-2*Py()),
					   Px()*exp(-2*Py())) );


    typedef bases<Lagrange<OrderChamp/*+3*/,Scalar,PointSetFekete> > basis_image_type;
    typedef FunctionSpace<mesh_type, basis_image_type> space_image_type;
    typedef boost::shared_ptr<space_image_type> space_image_ptrtype;
    typedef typename space_image_type::element_type element_image_type;

    typedef bases<Lagrange<OrderChamp/*+3*/,Vectorial,PointSetFekete> > basis_image_v_type;
    typedef FunctionSpace<mesh_type, basis_image_v_type> space_image_v_type;
    typedef boost::shared_ptr<space_image_v_type> space_image_v_ptrtype;
    typedef typename space_image_v_type::element_type element_image_v_type;


    space_image_ptrtype Yh = space_image_type::New( mesh2 );
    element_image_type u2( Yh, "u2" );

    space_image_v_ptrtype Yhv = space_image_v_type::New( mesh2 );
    element_image_v_type u2v( Yhv, "u2v" );

    typedef Backend<double> backend_type;
    typedef boost::shared_ptr<backend_type> backend_ptrtype;
    backend_ptrtype backend( Backend<double>::build( BACKEND_PETSC ) );

    //std::cout << "\nXh ndof = " << Xh->nDof();
    //std::cout << "\nYh ndof = " << Yh->nDof() <<"\n";
    /*_________________________________________________*/

    mesh->tool_localization()->updateForUse();

    boost::timer __t;
    __t.restart();

    //OperatorInterpolation<space_type, space_image_type> I( Xh, Yh, backend );
    //I.apply(u,u2);
    OperatorInterpolation<space_v_type, space_image_v_type> I( Xhv, Yhv, backend );
    I.apply(uv,u2v);
    //u2v=I(uv);

    double __time = __t.elapsed();
    std::cout<<"\n Time : " << __time <<"\n";
    //double err = integrate( elements(mesh), (idv(u)-idv(u2))*(idv(u)-idv(u2)) ).evaluate()(0,0);
    double err = integrate( elements(mesh), trans(idv(uv)-idv(u2v))*(idv(uv)-idv(u2v)) ).evaluate()(0,0);
    std::cout <<" \nerr : "<< std::sqrt(err)<<"\n";

    //dm.nGlobalElements(), dm.nMyElements()
    //std::cout <<"\n"<< u.map().nGlobalElements()<<" " << u.map().nMyElements()<<"\n";
    /*_________________________________________________*/

    //exportResults<mesh_type,element_type>( test_app, u, "oli" );
    //exportResults<mesh_type,element_type>( test_app, u2, "ole" );


  }


  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////

  /*_________________________________________________*
   * Options and About
   *_________________________________________________*/

  inline
  po::options_description
  makeOptions()
  {
    po::options_description desc_options("test_interp_twomesh options");
    desc_options.add_options()
      ("hsize1", po::value<double>()->default_value( 1 ), "mesh size1")
      ("hsize2", po::value<double>()->default_value( 0.5 ), "mesh size2")
      ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh(2.0 or 2.1)")
      ;
    return desc_options.add( Life::life_options() );
  }

  inline
  AboutData
  makeAbout()
  {
    AboutData about( "OP_Interp" ,
		     "OP_Interp" ,
		     "0.1",
		     " OP ",
		     Life::AboutData::License_GPL,
		     "Copyright (c) 2010 Universite Joseph Fourier");
    about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
    return about;
  }
}

int main(int argc,char** argv)
{

  Application_ptrtype test_app(new Application_type(argc,argv, opinterp::makeAbout(), opinterp::makeOptions() ));

  test_app->changeRepository( boost::format( "/opinterp/%1%/" )
			      % test_app->about().appName()
			      );

  const uint nDim = 2;
  const uint nOrderChamp = 5;
  const uint nOrderGeo = 1;

  opinterp::run<nDim,nOrderChamp,nOrderGeo>(test_app);

  return 0;

}
